Trapped fermion mixtures with unequal masses: a Bogoliubov-de Gennes approach 
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We use the Bogoliubov-de Gennes formalism to analyze the ground state phases of harmonically 
trapped two-species fermion mixtures with unequal masses. In the weakly attracting limit and 
around unitarity, we find that the superfluid order parameter is spatially modulated around the 
trap center, and that its global maximum occurs at a finite distance away from the trap center 
where the mixture is locally unpolarized. As the attraction strength increases towards the molecular 
limit, the spatial modulations gradually disappear while the Bardeen-Cooper-Schrieffer (BCS) type 
nonmodulated superfluid region expands until the entire mixture becomes locally unpolarized. 
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The many-body physics of fermion mixtures with mis- 
matched Fermi surfaces has been a longstanding problem 
for many researchers ranging from the condensed and nu- 
I clear matter to the high energy and astrophysics commu- 
nitites [ij. While there are various theoretical propos- 
als for the ground state of such systems including the 
Fulde-Ferrell-Larkin-Ovchinnikov (FFLO), Sarma and 
breached pair superfluid phases, strong experimental ev- 
idence for their observation is still lacking. Following the 
recent experiments on two-component ^Li mixtures with 
unequal populations @, Q j ^^"^ wave of theoretical in- 
terest in this problem has been sparked in the atomic and 
molecular physics communities. In these experiments the 
phase diagram of trapped (finite) systems has been stud- 
ied as a function of population difference, temperature 
and two-body scattering length, showing superfluid and 
normal phases and a phase separation between them Q ■ 

Motivated by these experiments, phase diagrams of 
harmonically trapped mixtures with unequal populations 
have been extensively analyzed in both three- and one- 
dimensional systems. At the mean-field level in three di- 
mensions, while fully quantum mechanical Bogoliubov-de 
Gennes (BdG) calculations provide some evidence for the 
FFLO type spatially modulated superfluid phase, such a 
phase is completely absent in calculations based on the 
semi-classical local density approximation (LDA) 0, Q . 
Therefore, it is still an open question whether these spa- 
tial modulations are related to the FFLO superfluidity 
or are simply finite size effects. However, in exactly 
tractable one dimensional systems, FFLO structure of 
the superfluid phase have been identified in trapped as 
well as infinite systems 

0, i, S [13, Eli- These works 
arguably suggest that the ground state of polarized mix- 
tures is also an FFLO type superfluid in three dimensions 
along with the earlier BdG results Q . 

Two-species fermion mixtures with unequal masses of- 
fer a very natural way of creating superfluidity with mis- 
matched Fermi surfaces, and there have been increas- 



tems. For instance, ^Li-^'^K mixtures have recently been 
trapped and interspecies Feshbach resonances have been 
identified [2^ . [23| . opening a new frontier in ultracold 
atom research to study exotic many-body phenomena. In 
this manuscript, we go beyond the LDA method 15, 1^, 
and use the BdG formalism to analyze the ground state 
phases of harmonically trapped ^Li-'^'^K mixtures. Our 
main results are as follows. In the weakly attracting 
limit and around unitarity, we find that the superfluid 
order parameter is spatially modulated around the trap 
center, and that its global maximum occurs at a flnite 
distance away from the trap center where the mixture 
is locally unpolarized. As the attraction strength in- 
creases towards the molecular limit, the spatial modu- 
lations gradually disappear while the Bardeen-Cooper- 
Schrieffer (BCS) type nonmodulated superfluid region 
expands until the entire mixture becomes locally unpo- 
larized as shown in Fig. [T] 




ing theoretical [12 . 
and experimental 22l I 




FIG. 1: (Color online) The local superfluid order parameter 
|A(r)| versus distance r is shown for a population balanced 
mixture of ''Li and ''"K atoms. The BCS type superfluidity 
first occurs at a finite distance away from the trap center 
where the mixture is locally unpolarized. As the attraction 
strength increases towards the molecular limit, it gradually 
expands until the entire mixture becomes locally unpolarized. 



We obtain these results by using the following Hamil- 
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tonian (in units of 7i = fc^ = 1) 



H ^ dr 



V4(r)/C,(r)V.(r)-5*^(r)*(r) 



, (1) 



to describe two-component fermion mixtures with zero- 
ranged attractive {g > 0) interactions, where V'J(r) and 
Tpa (r) field operators create and annihilate a pseudo-spin 
a fermion at position r. Here, we have introduced the 
operators ICa{r) = -S/^/{2M^) - /i^(r) and *(r) = 
V'j,(r)V'T(r)j where Ma is the mass, fia{r) — f^a — Vair) 
is the local chemical potential, fi^ is the global chem- 
ical potential and K-(r) = is the trapping 
potential of the a fermions. In the mean-field approxi- 
mation for the superfluid phase, this Hamiltonian reduces 
to Hmf = /drE,V'i(r)/C,(r)V.(r)-A(r)V|(r)Vl(r)- 
A*(r)-0|(r)7/'|(r) -I- |A(r)p/(7], where the self-consistency 
field A(r) = g{4>^{v)tpi{r)) is the local superfluid order 
parameter and (...) is a thermal average. 

Hmf can be diagonalized via the Bogoliubov- 
Valatin transformation '(/'^(r) — J2r]\-'^v,'^i^)^v-o- ^ 
^a-v^^ai^)'yl.~a-]^ where u^_cr(r) and Vri.a{j^) are the wave- 
functions and 7^ ^ and 7,,^cr are the operators correspond- 
ing to the creation and annihilation of a pseudo-spin a 
quasiparticle, and S| = +1 and S| = —1. This leads to 
the BdG equations 



/CT(r) A(r) 
A*(r) -/Ct(r) 



<y37),cr(r) = SCTe,,,CT<y5l,,cr(r), (2) 



where e^^^- > are the eigenvalues and ip^^^ij) are the 
eigenf unctions given by ipfj^^^r) — [u* |(r), w* ^(r)] and 
(fjj^l^r) — [v,j^i^{r),—Urj_i{r)]. Since the BdG equations 
are invariant under the transformation z)^^|(r) — > u* |(r), 
Ur^^l{r) —v* ^{r) and e^^| ~^n,h sufficient to 
solve only for w,,(r) = u^_-f(r), w,,(r) = w,,,x(r) and e,, = 
e^j as long as we keep all of the solutions with positive 
and negative eigenvalues. 

We assume A(r) = -g^^Ur^{Y)v'^{v)f{er^) is real 
which is sufficient to describe both nonmodulated and 
spatially modulated superfluid phases. Here, f{x) = 
l/[exp(x/T) -I- 1] is the Fermi function and T is the 
temperature. This equation has an ultraviolet di- 
vergence due to zero-ranged interactions, and it can 
be regularized by relating g to the two-body scatter- 
ing length ap of the t and J, fermions [2J]. This 
leads to 1/5 = -1/^/(4™^) + Mrkdr) / {2tt'^) - 
[M,kF{r)/{^n^)]\og{%{r) + kF{r)]/[k,{r) ~ fcf^(r)]}, 
where Mr = 2M^^MJ{M^ + M^) is twice the reduced 
mass of the t and | fermions, kp{r) / {2Mr) = [^t(^) + 



where na{r) = {'ipl{r)tpcr{j^)) is the local density of a 
fermions. This leads to ri|(r) — X^ij ^^'^ 

"i(r)=E,k.(r)lV(-e.). 

It is very natural to expand u,,(r) and w,,(r) in 
the complete basis of the harmonic trapping potential 
eigenf unctions, which are given by the Schrodinger's 
equation JCa{r)(f>^j^^{r) = ^^^^^^(r), where f = 
uja{2n + i + 3/2) — /icr is the eigenvalue and (j)'^ ^ m(r) = 
R'^j{r)Ye^jn{Sr, Vr) is the eigenfunction. Here, n is the 
radial quantum number, and £ and m are the orbital an- 
gular momentum and its projection, respectively. The 
angular part Yg^mif^r, 'i^r) is a spherical harmonic and 
the radial part is RZ^f,{r) = y2(Af^tj^)3/4[n!/(n + £ + 

l/2)!]i/2e-^'/2^^^1+'/'(r2), where = VM^r is di- 
mensionless and L^x) is an associated Laguerre poly- 
nomial. Since £ and m are good quantum numbers 
(77 = {£, m,7}), this expansion leads to u^,m,7(r) = 

E„C£,7,„(/)l^^^^„(r) and ^^.^.^(r) = E„ c'f,7,»<Ai,f,m('^)- 

The spherical symmetry of the Hamiltonian simpli- 
fies the numerical calculations considerably such that the 
BdG equations reduce to a 2{ne -I- 1) x 2{ne + 1) matrix 
eigenvalue problem for a given £ state 



E 



K 

a: 



-K 




) ,(3) 



where rii — {ric ~ £)/2 is the maximal radial quantum 
number and ric is the radial quantum number cutoff to 
be specified below. Here, the diagonal matrix element 
is K^'" = gSn,n' and the off-diagonal matrix element 

is A;*^"' = / r^drA{r)Rl i;{r)Rl,/r), where S,,j is the 
Kronecker delta. Furthermore, this procedure reduces 
the order parameter equation to 



A(r) = -g 



2£+l 



RLnir)Rl^^^,{r)f{ei,,), (4) 



and the local density equations to 



^ 47r 



ir)Ri,,n'i^)fi^in)^ (5) 
RLnir)Ri^,n'ir)fi-^i,ym 



where we introduced R^ ^ n{r) — Cf^-y.„Ji,„ 



i-r.nir) = Q,7,n-RL(^) and 
R^^ Jr) = d£,7,„i?i^f (r). Notice that the {2£ + 1) fac- 
tors in Eqs. (|4|), ^ and fG]) are due to the degener- 
acy of each £ state. Furthermore, the total number 
equations become TV^ = J2e,j,ni'^^ + ^)^'e,j,nfi'^in) and 
^ll{r)]/2 and fc2(r)/(2M,) = e, + [^i^{r) + ^i^{r)]/2. Here, = E« 7 n(2^ + l)'^L.„/(-Q,7)- These equations gen- 



Cc is the energy cutoff to be specified below, and our 
results depend weakly on the particular value of Cc pro- 
vided that it is chosen sufficiently high. Furthermore, 
the order parameter equation has to be solved self- 
consistently with the number equations — J drua (r) 



eralize the BdG formalism developed in Ref. [25| to the 
case with unequal masses, unequal chemical potentials 
and/or unequal trapping potentials. Having discussed 
the BdG formalism, next we analyze the ground state 
(T = 0) phases. 
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First, we analyze the noninteracting {g = or ap ^ 
0~) case. In this case, the discrete energy spectrum can 
be written as ^ = ujainp + 3/2) where Up = 2n + £ is 
the principal quantum number. Therefore, for a given Up, 
the orbital angular momentum £ ranges from 0, 2, Up 
when Up is even, and it ranges from 1,3, ...,np when 
is odd. Since the single pseudo-spin degeneracy = 
X;"''(2£+l) of each np level is D^^ = (rip + l)(np + 2)/2, 
we can introduce the Fermi level np.a that corresponds 
to the maximal value of the occupied Up states at T = 0. 
The condition = J2n''=o leads to — {nF,a + 
l)('^F,cr + 2)(nKcr+3)/6, and the energy eigenvalue ep.a- = 
ujfj{np fj + 3/2) that corresponds to the Up — np^ state is 
the Fermi energy of the a fermions. For sufhciently large 
np we notice that ep a and No- have a simple relation 
ep',^^u^{6N^)^/^. 

At T = 0, we can approximately calculate the po- 
sition r* where the local polarization density p{r) = 
|n|(r) — n|(r)| becomes zero p(?'*) = 0. Using 
LDA, we find na-{r) = kp ^{r)/ {Gtt'^), where kp^^ir) = 
Afg-Wo- (^"1^ ~ r^)^/^ is the local Fermi momentum and 
rp,a = {4:SN^y/^y y/MaU^ is the Thomas-Fermi radius of 
the cr fermions. This leads to = rp^i[M^uJi / {M^ui^ + 
M|W|)]^/^, which is an important length scale because 
the formation of BCS type Cooper pairs is most favored 
in the momentum space regions when the Fermi sur- 
faces of t and I fermions have minimal mismatch, i.e. 
kp.^iir^,) = kp,i{r^,). Therefore, when g 0+, the nonin- 
teracting mixture first becomes locally unstable against 
the BCS type superfluidity at r*, as can be seen in our 
numerical calculations which is discussed next. 




FIG. 2: (Color online) The local density n^(r) of ''Li (with 
points) and ''"K (with lines and points) atoms versus distance 
r is shown for a population balanced mixture. The inset shows 
the close-up densities near the trap edge. 

For this purpose, we solve the BdG equations [31 Q, 
([5]) and ([6]) self-consistently as a function of the dimen- 
sionless parameter l/(kpap) where kp is, specified be- 
low. In particular, we consider an equal population mix- 
ture of '^Li and '*°K atoms with 7V| = iV| and M| 
O.lSMj, and assume that both ^Li and ^"K atoms are 



trapped with equal trapping potentials V|(r) = Vi{r) 
such that M-\ui^ = ^1^1- We introduce a 'reduced' 
trapping frequency ojr via MrUj^ — M^uj^, and define 
an energy scale ep and two length scales {kp,rp) via 
ep = LOrlnp -t-3/2) = kj./{2Mr) = MrUj'^rj,/2. Notice 
that np = np^i = np^i when iV| = N^. Similarly, we 
define via Cc = a;,- (nc -I- 3/2). In our numerical calcula- 
tions, we choose np = 15 and ric — 65, which correspond 
to a total of = 1632 fermions and ec ~ 4:ep, respec- 
tively. Here, it is important to emphasize that we do not 
expect any qualitative change in our results with higher 
values of ni? and/or Uc, except for minor quantitative 
variations. 

In Figs. [Tl [5] and O we show the evolution of A(r), 
na-{r) and p{r) as a function of \/{kpap). For a weakly 
attracting mixture, while ^"^K atoms are in excess around 
the trap center and ^Li atoms are in excess close to the 
trap edge, they have similar densities only around r =^ 
K, 0.85rF.| w 0.61rF. Therefore, when \/{kpap) < 
—0.5, we find that A(r) is spatially modulated around 
r = 0, and that its global maximum occurs at r = 
where the mixture is locally unpolarized. The modula- 
tion period T w Q.Qrp of A(r) is approximately given by 
T ~ 2t: /\kp,'^ — kp^i \ w O.hrp, where kp^^r = \/2Mo7v- As 
g increases towards unitarity \/{kpap) — 0, we find that 
the amplitude of the modulations dramatically increases 
around r = 0. This is because both n|(r) and ni{r) are 
highest at r = 0, which effectively leads to stronger lo- 
cal interactions there since \/[kp{r)ap] increases with in- 
creasing density when ap <Q. Therefore, the maximum 
of A(r) eventually occurs at r = when the effective 
local interactions become sufficiently strong. 

These spatial modulations of A(r) have dramatic ef- 
fects on the local density of fermions causing pronounced 
modulations in n„{r) and p{r) close to r = as shown in 
Figs. [2] and [3l Since these modulations are significantly 
large, we hope that they can be observed in the future 
experiments. Further increasing g towards the molecular 
limit l/{kpap) ^ 0.5, we find that the spatial modula- 
tions gradually disappear and the BCS type nonmodu- 
lated superfluid region expands until the entire mixture 
becomes locally unpolarized. This is expected because 
the Fermi surfaces disappear in this limit, and therefore 
formation of the molecules does not require matching of 
the Fermi surfaces. 

We remark in passing that the BdG equations do not 
necessarily have a unique solution, and depending on the 
initial values of fXa- and A(r) that are used in the iterative 
approach, they often yield multiple solutions for a given 
set of parameters. In this manuscript, we show only the 
physical solutions which have lowest energy. Compared 
to the physical solutions, the unphysical ones have con- 
siderable qualitative variations around r = 0, i.e. both 
A(r) and n^{r) have more modulations. However, both 
the physical and the unphysical solutions have very sim- 
ilar qualitative structure around r = r*, where the BCS 
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type nonmodulated superfluidity occurs. 
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FIG. 3; (Color online) The local polarization density p(r) — 
|n|(r) — n|(r)| versus distance r is shown for a population 
balanced mixture of ^Li and ''"K atoms. While ®Li atoms are 
in excess near the trap edge, atoms are in excess close to 
the center. 

We emphasize that our results are based on the fully 
quantum mechanical BdG formalism, and that they are 
significantly different from the earlier results that are 
based on the semi-classical LDA method [ll, [l^. For 
instance, these LDA calculations suggest a sharp phase 
separation in the weakly attracting limit and around 
unitarity such that the BCS type superfluid phase is 
sandwiched between locally polarized normal fermions. 
Therefore, both A(r) and ncr{r) have unphysical discon- 
tinuities at the normal-BCS superfluid-normal interfaces 
which indicates breakdown of the LDA. This is because 
the LDA method excludes the possibility of a spatially 
modulated superfluid phase, which is one of the possible 
candidates for the ground state. A similar discrepancy 
between the BdG and the LDA method was previously 
discussed in the context of one-species fermion mixtures 
with unequal populations 0, [^. Thus, we conclude quite 
generally that the LDA method is insufficient to describe 
trapped fermion mixtures with mismatched Fermi sur- 
faces, and that it should be used with caution. 

Furthermore, our results for population balanced two- 
species fermion mixtures are qualitatively different from 
the recent works on one-species fermion mixtures with 
unequal populations 0, 0]. In the latter case, the BCS 
type superfluid phase occurs around r = 0, and A(r) 
is spatially modulated towards the trap edge where the 
mixture is locally polarized. Since both n-\'{r) and ni{r) 
are very low near the trap edge and the modulations have 
very small amplitudes, it may not be possible to observe 
them at experimentally attainable temperatures. How- 
ever, in our case, spatial modulations with large ampli- 
tudes occur around r = where both n| (r) and (r) are 
very high. These make two-species fermion mixtures very 
good candidates for the observation of spatially modu- 
lated superfluid phases in atomic systems. 

For instance, spatial modulations of A{r) can be ob- 



served by using the recently developed technique of spa- 
tially resolved radio- frequency spectroscopy 26]. This 
technique can be used to locate the nodes of A(r), since 
the local quasiparticle excitation spectrum becomes gap- 
less at the position of the nodes. In addition, spatial 
modulations of ^^(r) can be observed by using phase- 
contrast imaging of ^Li and ''"K populations. In fact, 
both of these techniques have recently been used with 
great success to characterize the superfluid and the nor- 
mal phases of one-species fermion mixtures with unequal 
populations [lH. 

In conclusion, we analyzed the ground state phases of 
harmonically trapped ^Li-^°K mixtures with equal pop- 
ulations. In the weakly interacting limit and around uni- 
tarity, we found that the superfluid order parameter is 
spatially modulated around the trap center. Further- 
more, we showed that the BCS type superfluidity first oc- 
curs at a finite distance away from the trap center where 
the mixture is locally unpolarized, and then it gradually 
expands as the attraction strength increases towards the 
molecular limit until the entire mixture becomes locally 
unpolarized. Since the spatial modulations with large 
amplitudes survive at unitarity, two-species fermion mix- 
tures offer a unique opportunity for their observation. 
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